Quick Mapping Lightning Talk

Setup

library(dplyr)
library(sf)
library(readr)
library(mapview)
library(Rspatialworkshop)
library(ggplot2)
library(leaflet)
library(httr)

Get flat data with coordinates

Use read_csv fromreadr package to load a .csv file I keep as part of a spatial data package that you can install with:

library(devtools)
install_github("mhweber/Rspatialworkshop")
library(Rspatialworkshop)

Then we:

1. Read in the csv
2. Suppress showing column types
3. Select a subset of numerous attributes
4. Filter a small portion of the 2700 gages in the dataset by filtering only gages with drainage areas > 500 sq miles

gages <- read_csv(system.file("extdata/Gages_flowdata.csv", package = "Rspatialworkshop"),show_col_types = FALSE) %>% dplyr::select(ID=SOURCE_FEA,LON_SITE,LAT_SITE, MeanFlowCFS=AVE, DrainAreaSqMiles=DA_SQ_MILE) %>% 
  dplyr::filter(DrainAreaSqMiles > 500)
glimpse(gages)

Rows: 539
Columns: 5
$ ID               <dbl> 14361500, 14378000, 14335000, 1…
$ LON_SITE         <dbl> -123.3178, -123.8123, -122.5956…
$ LAT_SITE         <dbl> 42.43040, 42.37900, 42.69985, 4…
$ MeanFlowCFS      <dbl> 3402.544, 2317.657, 1796.864, 7…
$ DrainAreaSqMiles <dbl> 2459, 665, 650, 698, 1610, 670,…

Make our data spatial

• We use st_as_sf from the sf package to promote the data to a simple feature collection by:
  • passing the longitude and lattitude fields from our data to use
  • providing a coordinate reference system - we specify this using an epsg code
  • Our epsg code is for unprojected NAD83
  • We can find coordinate reference systems (CRS) easily from spatialreference.org

gages <- st_as_sf(gages, coords = c("LON_SITE", "LAT_SITE"), crs = 4269, remove = FALSE)
ggplot() + geom_sf(data=gages)

A simple interactive map using mapview

mapview(gages)

Customize mapview

• Change the background
• Scale the symbol size by flow rate

m <- mapview(gages, map.types = c("Esri.WorldShadedRelief", "OpenStreetMap"), color = "grey40",cex = "MeanFlowCFS")
m

Add Web Map Service (WMS) Layers!

Web Map Services provide pre-rendered map tiles at different scales and are useful as background map layers. Here is my original Stack Overflow question from when I first played with adding WMS layers with mapview and leaflet

NHD WMS from the National Map

Here we display a background of National Hydrography Dataset stream lines to go with our gages without having to load local files.

This is a super handy way to share your data in an interactive map mashed up with relevant background layers in one file in R, without having to share spatial files.

m <- mapview(gages, cex = "MeanFlowCFS")
# Define the NHD WMS service
wms_nhd <- "https://hydro.nationalmap.gov:443/arcgis/services/nhd/MapServer/WmsServer?"
m@map <- m@map %>% 
  addWMSTiles(group = 'USGS HydroCache',
              wms_nhd,layers  = 6,
              options = WMSTileOptions(format = "image/png", transparent = TRUE),
              attribution = "")  %>% mapview:::mapViewLayersControl(names = "NHD")

m

Other WMS Services

We can see other services available on the National Map here - let’s try transportation

# Define the WMS service
m <- mapview(gages, cex = "MeanFlowCFS")
wms_trails <- "https://carto.nationalmap.gov:443/arcgis/services/transportation/MapServer/WmsServer?"
m@map <- m@map %>% 
  addWMSTiles(group = 'Trails',
              wms_trails,layers  = 3,
              options = WMSTileOptions(format = "image/png", transparent = TRUE),
              attribution = "")  %>% 
  addWMSTiles(group = 'Trail Labels',
              wms_trails,layers  = 15,
              options = WMSTileOptions(format = "image/png", transparent = TRUE),
              attribution = "") %>%  mapview:::mapViewLayersControl(names = c("Trails","Trail Names"))
                                  
m

Querying and making spatial Web Feature Services (WFS)

Just to show scope of different things we can do linking to REST Services as well - here are links to a couple examples- spatial wfs services and accessing arcgis rest services

And we can also use REST services for the same WMS layers we displayed in maps above! See National Map Services page

We can get a listing of all ESRI REST services here

Query water bodies REST feature

url <- parse_url("https://services.arcgis.com/P3ePLMYs2RVChkJx/arcgis/rest/services")
url$path <- paste(url$path, "USA_Water_Bodies/FeatureServer/0/query", sep = "/")
url$query <- list(where = "STATE = 'OR'",
                  outFields = "*",
                  returnGeometry = "true",
                  f = "geojson")
request <- build_url(url)


wb <- read_sf(request)

mapview(gages) + mapview(wb, col.regions='light blue')

Query trails we used earlier

This time use the REST service and a query for ‘Pacific Crest National Scenic Trail’ - bonus add a default topography brackground with it

url <- parse_url("https://carto.nationalmap.gov/arcgis/rest/services")
url$path <- paste(url$path, "transportation/MapServer/37/query", sep = "/")
url$query <- list(where = "Name = 'Pacific Crest National Scenic Trail'",
                  outFields = "*",
                  returnGeometry = "true",
                  f = "json")
request <- build_url(url)

PCT <- read_sf(request)

mapview(PCT, color='dark green',map.types = "OpenTopoMap")